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We investigate, mostly analytically, the nature of the Fulde-Ferrell- Larkin-Ovchinnikov super- 
fluid phases for an isotropic two-dimensional superconductor, which is relevant for the experimental 
observation of these phases. We show that there is, in the low temperature range, a cascade of 
transitions toward order parameters with ever increasing complexity. 
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The transition to Fulde-Ferrel-Larkin-Ovchinnikov 
(FFLO) superfluid phases [1,2] has given rise to much 
work, but a full theoretical understanding of these inho- 
mogeneous phases has not yet been reached. This is of 
high practical importance since experiment relies heavily 
on theory when it needs to test if these phases have actu- 
ally been observed. Indeed beside their basic theoretical 
interest which, in addition to superconductivity and its 
coexistence with magnetism [3], ranges from ultracold 
gases [4] to quark matter as it might occur in the core of 
neutron stars [5,6], FFLO phases are also of very high ap- 
plied interest, because they correspond to the supercon- 
ducting phases which should appear in superconductors 
with very high critical fields in the vicinity of the tran- 
sition to the normal state. A main strategy to observe 
these transitions in superconductors is to get rid of the or- 
bital currents, responsible for the standard disappearance 
of the superconducting phase at low critical fields, and 
reach the regime of Pauli paramagnetic limitation. This 
can be done in clean quasi two-dimensional systems such 
as organic compounds made of widely separated conduct- 
ing planes, or in high T c compounds. In these cases hop- 
ping between planes is very strongly inhibited. So when 
a strong magnetic field is applied parallel to the planes, 
the resulting orbital currents perpendicular to the planes 
are very weak and there is basically no orbital pair break- 
ing effect which leaves open the possibility of transitions 
to FFLO phases at much higher fields. Nevertheless the 
small intcrplanc coupling is necessary to suppress phase 
fluctuations and produce phase which can be properly 
described by a mean field theory. Indeed observations of 
this transition have been claimed quite recently in such 
organic compounds as k-(BEDT- TTF) 2 Cu(NCS) 2 [7,8] 
or A- (BETS- TTF) 2 GaCl 4 or FeCl 4 [9,10]. 

Actually when the field is not exactly parallel to the 
planes there are small orbital currents and associated vor- 
tices appear. The physics of the evolution from this vor- 
tex state to the pure FFLO phase is quite remarkable, 
both in the linear [11-13,3] and in the non linear regime 
[14] since it involves a series of vortex states phases corre- 
sponding to Landau levels with increasing quantum num- 
ber. Compared to this complex structure the limiting 2D 



FFLO phase sounds simple. Indeed the FFLO transition 
in 2D systems is believed to be second order and in par- 
ticular Burkhardt and Rainer [15] have studied in details 
the transition to a planar phase, where the order parame- 
ter A(r) is a simple cos(q.r) at the transition. This phase 
has been found by Larkin and Ovchinnikov [2] to be the 
best one in 3D at T = for a second order phase transi- 
tion. And in 3D it is also found to be the preferred one 
in the vicinity of the tricritical point and below [16-18], 
although in this case the transition turns out to be first 
order (except at very low temperature). However it is 
not clear that this simple order parameter is always the 
best one since, as first explored by Larkin and Ovchin- 
nikov, it is in competition with any superposition of plane 
waves, provided that their wavevectors have all the same 
modulus. 

In this paper we explore the low temperature range 
in 2D for the pure FFLO phase transition in the sim- 
ple isotropic case and show that surprisingly there is a 
cascade of the second order transitions toward order pa- 
rameters with ever increasing complexity. In principle 
there is an infinite number of phases. Therefore the pure 
FFLO phase display a structure just as rich as the vortex 
phases that we mentionned above, and the interplay be- 
tween both promises to be of exceptional interest. A first 
step in this direction has been made recently by Shima- 
hara [19] who, instead of simple cosine, found a tran- 
sition toward a superposition of up to six plane waves 
when looking for more complex 'cristalline' structures. 
We find agreement with his results. However we will see 
that cristallinity does not give the proper understanding 
of these phases. Instead we show that, when the temper- 
ature is lowered toward T = 0, the cascade of transitions 
toward order parameters with an ever increasing number 
of plane waves is due to the singular nature of the T = 
limit, which is already appearent in the T = wavevector 
dependence of the second order term in the free energy 
expansion in powers of the order parameter. 

Our study follows the general framework set up by 
Larkin and Ovchinnikov [2] and used by Shimahara [19] 
in the 2D case. The free energy difference = fl s — f2 n 
between the superconducting and the normal state is ex- 
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pandcd to fourth order in powers of the order parame- 
ter A(r) = Aq. exp(iqj.r). The second order contri- 
bution to is N J^q ^2(9, M, T)|A q | 2 , where N is the 
single spin density of states at the Fermi surface and 
p, = (/i-j- — /ij.)/2 is half the chemical potential differ- 
ence between the two spin populations. At the FFLO 
transition we have ^2(9, fiiT) = 0. Naturally the actual 
transition corresponds to the largest possible ft at fixed 
T. It is easily seen [2] that the corresponding wavevector 
q for the order parameter is obtained by minimizing: 
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with respect to q, where we have introduced the di- 
mensionless wavevector q — qvp/2p, and u) n = u) n — ifi 
with u n — irT(2n + 1) being the Matsubara frequency. 
This can be rewritten as an integral over real frequency. 
At T = one finds explicitely I(q,fi,T) = Reln(l + 
\J\ — q 2 ) — In 2, which is minimum for q — 1 , in agree- 
ment with Shimahara [20] and Burkhardt and Rainer 
[15]. We note that the location q = 1 of the minimum 
is a singular point for / since / = ln(g/2) for q > 1, so 
I is continuous while its derivative is discontinuous for 
q = 1. For T^O but low enough temperature, the con- 
dition giving the minimum can be expanded in powers 
of t 1 / 2 where t = T/fi is a reduced temperature. We 
have found that expanding q— 1 up to t 2 gives a result in 
excellent agreement with a straight numerical evaluation 
for t < 0.1. Here (and similarly below) we do not give 
the somewhat lengthy full expression [21], but we omit 
the last two terms and display only the lowest order cor- 
rection which is: 



1 = ^ 
2 2t 



(2) 



The second order term in Q fixes only the modulus, 
given by Eq. (2) , of the N wavevectors entering the plane 
wave expansion of the order parameter. This leaves open 
all the various possible order parameters arising from any 
combinations of these plane waves. The selected state 
will correspond to the lowest fourth order term, which is 
given [2] in general by: 
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^(2-,5 qi , q3 )|A q J 2 |A qj | 2 JK i , q3 ) (3) 



+(1 - <5 qs , q , - <5 qi ,_ qj )A qi A_ qi A q .Al q . J(a qiiqj ) 

where ctq ii q j is the angle between and qy. Just as for 
the second order term J(a) and J(a) can be expressed 
[21] as frequency integrals, expanded at low temperature 
t ^ 1 and evaluated for the optimal value of q (given by 
Eq.(2) to leading order). As for the optimal q we have 
carried this expansion up to second order in t 1 / 2 . If we 
restrict ourselves to leading order as in Eq.(2) we find for 
example for J (a): 



16tfl 2 J(a) = - 



+ cos( 



a cosh 2 [(1/4) ln§-/3 2 /2] 

r, exp(-^ 2 ) 

a/2)) Jo v 2 + P 2 



(4) 



where we have set (3 2 = (1 — cos(a/2))/i. The result up 
to second order in t 1 ! 2 is not much more complicated. A 
similar result [21] is found for J (a). 
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FIG. 1. J(a) for values 0.01 and 0.05 of the reduced tem- 
perature i = T/(qfl) = t/q (with essentially q ^ 1 in all this 
range). The lines are calculated using our low temperature 
expansion for J (a), while the points correspond to the di- 
rect numerical summation over Matsubara frequencies. The 
dashed line in the insert (which displays the same results with 
different scale) is the asymptotic behaviour for J(a) in the 
T -> limit Eq.(5). 

Plots of J(a) for various low temperatures can be seen 
in Fig. 1, either obtained from these low temperature ex- 
pressions or from exact calculations of J(a) obtained by 
direct numerical summation over Matsubara frequencies 
after analytical angular averaging. We see that, at this 
level of accuracy, both agree remarkably well up to rather 
high temperatures which shows that our low tempera- 
ture expansion is quite under control. Now Fig.l shows 
that J(a) has at low T a quite remarkably structured be- 
haviour which can be easily understood from some lim- 
iting cases. First let us take the limit T — > at fixed a. 
This implies (3 2 — > 00, the first term in Eq.(4) goes to 
zero and the easy integration leads to: 



J (a) 



4^ 2 sin 2 (a/2) 
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On the other hand if at fixed T we take the limit a —> 0, 
we have [3 2 — > 0. The limiting behaviour of the integral 
is it /(2(3) — 7T 1 / 2 , the dominant divergent contributions 
from the two terms in Eq.(4) cancel out giving: 
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which goes naturally to infinity for T — > 0. From these 
two limiting situations we can understand that, for most 
of the a range, J(a) is negative as it can be seen from 
Eq.(5) and it goes to large negative values when a gets 
very small. On the other hand for a — or very small 
J (a) is positive and very large, as results from Eq.(6) 
(surprisingly J(a) starts first to increase strongly before 
going down to very negative values). 

We have made the same kind of treatment [21] for 
J {a). The limit T — > with fixed a leads to p, 2 J(a) = 
— 7r/4sina while in the limit a — > at fixed T we find 
jl 2 J(a) = —1/4. These two cases make reasonable that 
J(a) is always negative, as we find. However, just as 
for J(a), one finds also a singular behaviour at small 
a. While for (3 w a/ {St) 1 / 2 » 1 and p?J(a) diverges as 
— 7r/4a, it goes to the finite value —1/4 for a = 0. Never- 
theless the divergent behaviour in a -1 is weaker than the 
one in a~ 2 found for J (a) and similarly J(0) » |J(0)| at 
low T. So J (a) will play the dominant role and at first 
we omit J (a) from our considerations. 

Let us now consider which is the most favorable or- 
der parameter. We will always find that the fourth order 
term is positive, so we have to make it as small as pos- 
sible in order to minimize the free energy. We note first 
that there is unevitably a strongly positive contribution, 
proportional to J(0), coming from the N terms = in 
the first sum of Eq.(3). However their unfavorable effect 
is relatively reduced when N increases since there are TV 2 
terms in this sum. Since we clearly have to avoid other 
positive contributions of this kind, we want to forbid the 
angle domain < a < ao where we can take ao as giving 
J(ao) = 0, which is very close to the (strongly negative) 
minimum of J (a) as it can be seen on Fig. 1. So the an- 
gle between any two different wavectors q, and q^ should 
be at least ao . On the other hand it is of interest to have 
angles close to ao since they give strongly negative con- 
tributions to O. Since it is better to have the number of 
wavevectors TV as large as possible, we come from symme- 
try reasons to the conclusion that we have to choose the 
wavevectors qi angularly equally spaced with the angle 
between two nearest wavevectors as close as possible to 
a . Now a — > when T — > 0. Therefore we have indeed 
a cascade of transitions corresponding to order param- 
eters with an ever increasing number N of wavevectors 
when T — > 0. This singular behaviour is clearly linked to 
the fact that it is not possible, in particular for J (a), to 
perform an expansion in the limit T — > 0, i.e. this limit 
is singular. 

More precisely it is quite a difficult problem to han- 
dle this minimization from scratch, and for example to 
prove that the minimum is not obtained for a com- 
pletely disordered situations. However symmetry makes 
such a result quite unlikely. On the other hand if we take 



for granted that the wavevectors are equally spaced an- 
gularly it is easy to prove that the weight |A q J are all 
equal, because this problem can be mapped on a tight 
binding Hamiltonian on a ring. Conversely if we assume 
that all these weights are equal the problem is equiva- 
lent to find the equilibrium position of atoms on a ring 
with repulsive short range interaction and attractive long 
range interaction. We expect the equilibrium to be a 
crystalline structure corresponding in our case to equally 
spaced wavevectors. 
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FIG. 2. G^ 1 as a function of temperature for various plane 
waves number N (for clarity only solutions up to N = 14 are 
presented). The arrows give the values of the successive tran- 
sitions. In the upper panel the dashed lines are the exact Mat- 
subara summations. The full lines in the upper panel, as well 
as all the lines in the lower panel, are from the low tempera- 
ture expansion for J (a) and J(a). In the insert the full dots 
are the exact temperature locations of the transitions together 
with the corresponding values of N (the empty dots just give 
N — 2). The full line gives the low temperature evaluation of 
these transitions from the solution of a simple transcendental 
equation, and the almost undistinguishable dashed line comes 
from the explicit first iteration solution of this equation. 
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In the low temperature limit (t -C 1), we can easily be 
more quantitative and write simple analytical answers for 
our problem. Indeed the result is dominated by the small 
a behaviour of J {a). In this range Eq.(4) becomes: 



H 2 J{a) ~ 



1 



16iacosh X 



(7) 



where X = x - (1/4) ln(7r/2i) with x = a 2 /IQt. The 
second term is the asymptotic form Eq.(5). The zero of 
J(a) we are looking for is found for large X and sat- 
isfies e x = (k/y/2)x^ 4 with k = (2n 3 /t 2 ) 1 ^, which is 
easily solved iteratively for large k. The leading order 
x = ln(fc/V2) gives a = {St ln[(7r 3 /2) 1/2 /i]) 1/2 . Thc 
first order iteration is x = ln(fc/\/2) + l/41n(ln(fe/\/2). 
Then the optimum value of AT is given by the integer 
value of 27r/ao, which conversely fixes the temperature 
at which the system switches from N — 2 to iV wavevec- 
tors. In the insert of Fig. 3, one sees that the first iteration 
gives results almost undistinguishable from the exact so- 
lution of the above transcendental equation for x. 

Exact results can be obtained by handling numerically 
the various steps. The minimization of the free energy 
with respect to |A q | gives fl/N = -fl 2 (q, fi,T)/G 2 (N) 
with: 



N-l 



NG 2 (N) = 2J(0) + 4 J(^n/N) 



(8) 
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where we have omitted the J terms for simplicity. One 
sees in Fig. 3 the results (including the J terms) for 
G^ 1 (AT), i.e. essentially the free energy for various values 
of N. The stable state has the highest (A7). The cas- 
cade of phase transitions is clearly seen in the lower panel, 
corresponding to the lower temperature range. The ar- 
rows show the successive transition temperatures. The 
results from straight numerical Matsubara summation 
can not be distinguished from those obtained from our 
low T analytical expression for J (a) and J {a), except for 
the higher temperature, as seen in the upper panel. The 
critical temperatures for switching from N — 2 to N are 
shown as dots in the insert of Fig. 3. It can be seen that 
our above asymptotic calculation for these temperatures, 
shown as full curves, are actually quite good. 

At low temperature we can make an asymptotic calcu- 
lation of G2(N). As seen on Fig. 1 J (a) switches rapidly 
above ag to its large angle asymptotic behaviour. More- 
over at low T the number of plane waves N is large and in 
Eq.(8) the sum is dominated by the small angles terms, 
which leads [21] to: 



J1 2 G 2 (N) ~ 



2Nt 



N 
~3~ 



(9) 



This expression corresponds in Fig. 3 to the rising part, 
on the low temperature side, of G 2 ~ 1 (A'). On the other 
hand the downturn is due to contributions from the pos- 
itive part of J (a) which can not be evaluated so easily. 



When we substitute in Eq.(9) the value N — 2-ir/ao for 
the optimum plane wave number we have: 



fi 2 G 2 (N) 



4irt 
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Since St/ot^ ~ l/ln(l/t) from our above evaluation, we 
find that G 2 (N) is always be positive in the low tempera- 
ture range, which means that the transition stays always 
second order. When we take into account the contribu- 
tion of the J terms in Eq. 3 that we had omitted so 
far, one gets an additional — In N contribution in Eq.(9), 
which arises only when N is even because for equally 
spaced wavevectors, the J terms contribute only in this 
case since they require opposite wavevectors. This con- 
tribution is just enough to systematically tilt the balance 
in favor of even N. In this respect the phase N = 3 found 
by Shimahara [19] appears as an exception. Otherwise 
our cascade of phase transitions contains only order pa- 
rameters with an even number of plane waves, and the 
order parameter is just a real sum of cosines, with still a 
degeneracy due to the freedom of choosing the phases in 
these cosines. 

We note that, although the transitions between the 
various FFLO phases and the normal phase are second 
order, we expect the transitions between different FFLO 
phases in the superfluid domain to be first order since 
one can not go continuously from a given order param- 
eter to the next one. Finally we have seen that these 
successive transitions are directly due to the singularity 
which occurs in 2 D at T = 0. This singularity itself 
arises because the two Fermi circles corresponding to op- 
posite spins come just in contact when one applies the 
shift corresponding to the wavevector q of the FFLO 
phase. Therefore this singularity is a general feature 
of 2D physics and we expect it to give rise to similar 
consequences in more realistic and more complex models 
describing actual physical systems. 

* Laboratoire associe au Centre National dc la 
Recherche Scicntifiquc ct aux Universites Paris 6 et Paris 
7. 
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